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Abstract 

A biological signal is transmitted by interactions between signaling molecules in the cell. To date, there have been extensive 
studies regarding signaling pathways using numerical simulation of kinetic equations that are based on equations of 
continuity and Fick's law. To obtain a mathematical formulation of cell signaling, we propose a stability kinetic model of cell 
biological signaling of a simple two-parameter model based on the kinetics of the diffusion-limiting step. In the present 
model, the signaling is regulated by the binding of a cofactor, such as ATP. Non-linearity of the kinetics is given by the 
diffusion fluctuation in the interaction between signaling molecules, which is different from previous works that 
hypothesized autocatalytic reactions. Numerical simulations showed the presence of a critical concentration of the cofactor 
beyond which the cell signaling molecule concentration is altered in a chaos-like oscillation with frequency, which is similar 
to a discontinuous phase transition in physics. Notably, we found that the frequency is given by the logarithm function of 
the difference of the outside cofactor concentration from the critical concentration. This implies that the outside alteration 
of the cofactor concentration is transformed into the oscillatory alteration of cell inner signaling. Further, mathematical 
stability kinetic analysis predicted a discontinuous dynamic phase transition in the critical state at which the cofactor 
concentration is equivalent to the critical concentration. In conclusion, the present model illustrates a unique feature of cell 
signaling, and the stability analysis may provide an analytical framework of the cell signaling system and a novel formulation 
of biological signaling. 
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Introduction 

Protein interaction is essential for cellular activities such as 
cytoskeleton formation [1,2], protein assembly [3], and cellular 
signaling [4]. The MAPK signal pathway is an example of a 
cellular signaling pathway that has been extensively studied [5-8] . 
A transient binding cofactor such as ATP/GTP, or phosphory- 
lation of amino acid residues, controls signaling molecule 
interactions and the subsequent modification of signaling mole- 
cules. Resulting reaction cascades operate to transmit cellular 
signals [6] . In cell signaling, oscillation has been reported in many 
studies, with circadian rhythms being a well-known example [9]. 
Calcium ion signaling oscillation is another well-known phenom- 
enon [10]. Mathematical models have been proposed to explain 
signaling kinetic behaviors based on a set of kinetic equations. 
Systems biology approaches have also been developed in recent 
years [1,2,11]. 

Systems biology can describe the kinetics of a signal pathway 
usingsimultaneousequations of a complex reaction network. On 
the other hand, there are other types of models consisting of many 
simultaneous reaction rate equations, including more than ten 
variables [12,13]; furthermore, systems biology models including 
more variables are also known [14]. Signaling networks frequently 
include non-linear reaction autocatalytic processes. To date, there 



have been many fascinating models of such autocatalytic reactions, 
enabling bifurcation and/or bi-stability in association with 
physical theory [5,15]. However, autocatalytic models or positive 
feedback are not necessarily applicable to all biological signaling 
pathways. 

In the current study, to understand biological signaling 
pathways, we used the following three novel perspectives (A)-(C) 
based on a non-linear and non-equilibrium kinetic model, which 
included only two concentrations of the signaling molecule to 
describe the biological signaling pathway. (A) Ari equation of the 

continuity of the chemical concentration of c, (i = 1 , n) 

including chemical reaction items can be described using diffusion 
coefficients D„ kinetic coefficients ki, and concentrations of 
individual compounds c, as follows: 

rlr " " 

^=AV 2 f ,+ ^>r,-+ £ k t c t c j+ --- (1.1) 

In the above formula, the diffusion rate is hypothesized to obey 
Fick's law. In general, the diffusion items and chemical reaction 
items are thus separately described. On the other hand, because 
the biological signaling pathway network, including protein 
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interactions, is limited by the slow diffusion rate of the signaling 
molecular proteins, kinetic coefficients generally depend on the 
diffusion coefficient. Therefore, diffusion items and reaction items 
cannot necessarily be separated without validation in the biological 
reaction. In the current model, we therefore described kinetic 
coefficients in a diffusion-coefficient-dependent manner. 

(B) A feedback process due to non-linear self-catalytic reactions 
was not assumed in the current model, but instead, interactions 
between signaling molecules in their diffusion was assumed to give 
non-linearity to the model. 

(C) A model system far from equilibrium due to a continuous 
supply of chemicals from the outside was hypothesized. The main 
issue is how minimal extracellular changes can be transformed into 
intracellular environmental changes. We aimed to evaluate the 
behavior of the model around the critical state by perturbation 
expansion using a minimal change of the supplied molecule 
concentration. By this mathematical evaluation, we aimed to 
illustrate the dynamic continuous oscillatory concentration change 
of signaling molecules from a static state. 

In the current study, we constructed a novel model and aimed 
to evaluate the general intrinsic properties underlying cellular 
signaling based on signaling molecule interaction kinetics. Previous 
systems biology models have not necessarily focused on the 
diffusion process of proteins. Given the non-linearity during 
diffusion, we assumed kinetic instability of the signaling molecule 
interaction, and the sensitivity of the cell signaling in response to 
the environmental change was evaluated. The model system 
consists of several steps as follows: (i) the signaling molecule 
achieves an interaction active state by reversibly binding a cofactor 
that provides the signaling molecule with interaction activity; (ii) 
the signaling molecule has the ability to hydrolyze the cofactor; (iii) 
the signaling molecule interaction activity becomes lower when 
binding a hydrolyzed inactive cofactor compared to the signaling 
molecule binding an active cofactor; (iv) the signaling molecule has 
the ability to exchange the inactive cofactor with an active one; (v) 
active cofactors are supplied continuously from the outside. Thus, 
we set the interaction activity to be self-limiting, causing dynamic 
instability of the signaling molecule interaction. In the present 
model, we assumed that a signaling protein diffuses relatively 
slowly in the cytoplasm, and the whole signal transduction is a 
diffusion-limited reaction (assumption A.l). Following the protein 
interaction, one of the signaling molecules is phosphorylated or is 
bound to a cofactor such as GTP or ATP. In the kinetics of the 
protein interaction, fluctuation analysis was not fully performed in 
spite of the greater fluctuation in concentration relative to the 
solution of small molecules. Furthermore, we systematically 
analyzed the roles of the fluctuation in cellular signaling. 

Materials and Methods 

Numerical simulation 

Numerical calculations were performed using Mathematica 8 
(Wolfram Research, Inc., Champaign, IL). 

Results 

Protein interaction kinetics 

The model scheme is shown in Figure 1 . There are two types of 
signaling molecule, an active cofactor-binding signaling molecule 
(X), and an inactive cofactor-binding signaling molecule (Z). An 
active cofactor is non-hydrolyzed, and the inactive cofactor is the 
hydrolyzed type. X has the higher interactive activity and Z has the 
lower interactive activity. First, X can associate with oligomeric 
enzyme complex R consisting of X and Z, which transforms the 



active form X into the inactive form Z by hydrolysis of the binding 
cofactor and is to be released as Z irreversibly: 



X + R *=t R *-k-\, kinetic coefficient) 



R^Z + R (k 2 ) (2.1) 

In the above formula, we assumed that signaling proteins diffuse 
relatively slowly in the cytoplasm and that the dissociation rate 
(relating to the irreversible orientation k-i) of an encountering 
pair X and Z is significantly slower than the hydrolysis rate of the 
active cofactor changing into inactive cofactor (relating to k 2 ) 
(assumption 2; A.2). On a simple consideration of the diffusion 
limited step, when the kinetic rate can be described according to 
Fick's law using diffusion coefficients D x and D z of X and Z, 
respectively: 

ki=4n(D x +D z )/2 (2.2) 

Next, Z recovers its interaction activity by exchange active 
cofactor P into inactive cofactor P\ returning to X (Figure 1). 

Z + P^X + P' (fc 3 ) (2.3) 

Signaling molecules have the potential to hydrolyze the cofactor 
by interacting with identical species: 

X+X^X+Z (fc,) (2.4) 



X + Z^2Z (k 5 ) (2.5) 

In (2.4) and (2.5), likewise in (1.2): 

k 4 =4nD x (2.6) 



k 5 =4n(D x +D z )/2 (2.7) 



Kinetic equations of interaction signaling molecules 

Here, kinetic equations were set according to the above simple 
reaction cascade. The equations consist of protein interactive items 
and an item of small-molecule-cofactor exchange. When the 
concentration of the protein is sufficiendy small, the dependency of 
the diffusion coefficient on the concentration is linear [11,13—18]. 
In comparison with the exchange kinetic rate of the cofactor (2.3), 
the rate of macromolecular protein interaction that depends on 
the diffusion step can be regarded as significantly smaller ((1.1), 
(2.4), & (2.5)) as a general (assumption 1; A.l). In this case, the 
whole reaction system can be regarded as diffusion-limited 
[4,19,20], and the diffusion rate is given according to Fick's law 
using a gradient of concentration. Thus, we get the protein 
interaction kinetics equations using the diffusion coefficient: 
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Figure 1. Scheme of signaling molecule interaction. Individual globules or oblongs represent signaling molecules X, Y, Z, and receptor ft. 
Kinetic coefficients are shown next to the arrows. Outside and inside signify the outside and the inside of the cell, respectively. 
doi:1 0.1 371 /journal.pone.01 0291 1 .g001 



X = - k\ RX + k^PZ - k 4 X 2 - k 5 XZ 



(2.8) 



k 5 =D 5 cc 



D x + D z 



(2.10) 



Here, the exchange rate of the cofactor is expressed by the item 
k 3 ' PZ. Further, 



Using (2.10), (2.8) and (2.9) are rewritten as follows: 



Z = k 2 R+k]PZ + k 4 X 2 + k 5 XZ 



(2.9) 



X= -D 1 RX+pZ-D 4 X 2 -D 5 XZ 



(2.11) 



Because kinetic coefficients depend on diffusion coefficients, we 



set: 



Z = k 2 R +pZ + D 4 X 2 + D 5 XZ 



(2.12) 



k\ =D\ gc 



Dx+Dr 



Here we set the oligomer concentration as constant because 
de novo asssembly is considered to be much slower than monomer 
interaction at the steady state (assumption 3; A.3): 



p = k 3 'P 



R = DiRX e -k 2 R = 0 



(2.13) 



k 4 =D 4 azD x 



Here, a small letter affixed signifies values at the steady state. 
This assumption is based upon on the steady state in the protein 
assembly [21-24]. 
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Setting the right hand sides of Eqs. (2.11) and (2.12) gives: 



Fluctuation of diffusion coefficient 

Subsequently, let us consider the fluctuation of participant 
proteins. We set: 

X = X e + x, Z = Z e +z (3.1) 



dD\=ax — bz (3.5) 

Here, an increase in X contributes to a decrease in diffusion 
coefficient Dj in thefluctuation item, ax («>0), because of the 
higher interaction activity that reducesdiffusion; in contrast, 
increased Z contributes to increases in the diffusion coefficient 
in the fluctuation item bz (b>0). 

Likewise in (4.5), 

dD 4 = cx — dz (c,d>0) 
dD 5 =ex-fz (e/>0) 



In actuality, receptor R interacts with other proteins, X and Z, 
in the course of diffusion (Figure 1). In actual signaling pathways, 
signaling molecules associate with other signaling molecules and 
phosphorylate them or are phosphorylated by them. The diffusion 
coefficients can be altered in proportion to the signaling molecule 
concentration. By using the Gibbs-Duhem expression, the 
diffusion coefficient D of one macromolecule in the solution can 
be generally written as [13]: 



D- 



k B T 

n 



(l-N A vc/M)(l+2AMc + ■■ 



(3.2) 



where T is the temperature of the solution, k B is the Boltzmann 
constant, and rj is the frictional coefficient of a given macromol- 
ecule in solution. A is the second virial coefficient, v is the partial 
specific volume of protein with molecular weight M, and N A is 
Avogadro's number. The small letter c denotes the concentration 
of the solute. 

Further, we hypothesized that D R> D x , and D z , the diffusion 
coefficients of R, X and Z, are given by extension of (3.2) to the 
mix solution of two macromolecules, X and Z: 

D X {X,Z) = 

^(1-^X-^Z)(1 + 2A X M X X + 2A Z M X X + ...) 
n x \ Mr M z J 

D Z {X,Z) = 

fcgTA N A v x ' x 

nz V M x 

D R (X,Z) = 



NrV Z ' , 
M z 



(3.3) 

Z ){\+2A x 'M x X+2A z 'M x X + ■ ■ ■) 



^L(l- 1 ^X- ^Iz)(l+2A X "M X X+2A Z "M X X+ ■ 



1r V 



My 



M 7 



where v x and v z are the partial specific volumes of X and Z with 
molecular weights M x and M z , respectively. A X A X ', A x ", A z , A z , 
and A z " are the second virial coefficients. In actuality, X and Z are 
the same molecules except with bound ATP or ADP. The 
fluctuation of the diffusion coefficient is given by: 



dD\ ccd 



D X +D R 



+ 



dx 

3 

dz 



Dx+Dr 
D X + D R 



(3.4) 



And therefore, we set: 



Using (3.5) and (3.6), Eqs. (2.5) and (2.7) give the fluctuation 
kinetic equations: 



x = - {R(Di - aX e ) + 2D A X e + D 5 Z e }x 
+ (Ra -D 4 + 2cX + eZ)x 2 
+ (p- bX e - D 5 X e - dX 2 -fX e Z e )z 
- (Ds + Rb~ eX e + f e Z e )xz -fX e z 2 



(3.7) 



z = (2D 4 X e + D 5 Z e - cX 2 - eX e Z e )x + (D 4 - 2cX - eZ)x 2 

+ (D 5 X e + dX 2 +fX e Z e -p)z (3.8) 
+ (D 5 + 2X e d - eX e +f e Z e )xz +fX e z 2 

Further, using matrix formulation, 



x 




X 






= L 




+ 


z 




z 





(3.9) 

' (Ra -D 4 + 2cX + eZ)x 2 - (D 5 + Rb - eX e +f e Z e )xz -fX e z 2 
(D 4 - 2cX - eZ)x 2 + (D 5 + 2X e d - eX e +f e Z e )xz +fX e z 2 

Here, 

(3.10) 

-R(Di -aX e )-2D 4 X e -D 5 Z e p-bX e -D 5 X e -dX 2 -fX e Z e 
2D 4 X e + D 5 Z e - cX 2 - eX e Z e D 5 X e + dX 2 + fX e Z e -p 



Calculus simulation of concentration oscillation 

The time-course of the signaling molecule concentrations was 
simulated via the substitution of appropriate numerical values into 
(3.9). A numerical calculation was performed over a sufficiently 
long period to evaluate the trend of signaling protein behavior. In 
the current simulation, the fluctuation coefficients a, b, c, d, e, and 
/, of Dj, T>4, and D$ are of the same order of magnitude, 10 2 
(100~856). The concentrations of X and Z at the steady state are 
given by Eq. (2.14). On the basis of A.3, Dj, the diffusion 
coefficient of the assembling rate of X to R, is significantly smaller 
than D 4 and D s , which are diffusion coefficients of the assembling 
rate between X and X, and X and Z. Applying the above 
conditions, the simulation results are shown in Figure 2. When p 
increases to the values that satisfy: 
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(a) (b) (c) (d) 




Figure 2. Time-course of the fluctuation of the signaling molecules displays a chaos-like oscillation. Diffusion of active cofactor binding 

signaling molecule (X) and of inactive cofactor binding signaling molecule (Z). The Appendix S1 presents the simulation parameters, with the notation 

of Eqs. (3.9). p is (a) 0.795, (b) 0.81, (c) 0.84, (d) 0.88, (e) 0.96, (f) 1.00, (g) 1.12, and (h) 1.16. The upper graph shows two parametric plots of X, and Z. 

Red, and blue lines in the lower graph represent the concentrations of X, and Z, respectively. The horizontal axis represents time (0 £ f £ 200) and the 

vertical axis represents the concentrations of X, and Z, respectively. When p exceeds 0.80, chaos-like oscillation is observed. Mathematica cord when 

p = 0.795 (a) is shown below. Below is the simulation program when p= 1.0253: 

D1=0.28 

k2 = 0 . 00034580 

a = 800 

b = 656 

c = 100 

d=100 

e = 100 

f = 100 

p = 1 . 0253 

D4 = 1 56 

D5 = 156 

R=1 

X = k2/Dl 

Z = (k2 (D1'2 R+ D4 k2) ) / (Dl (D1 p - D5 k2) ) 

NDSolve[ { Derivative[1][ x] [ t] = = - (R (D1 - a X) +2 X D4+ D5 Z) x[ t] + (R a - D4+2 c X + e Z) x[ t] 2+ (p - D5 X - b X - d X'2 - f X 
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Z) z[ t] - (D5+ Rb - e X + f Z) x[ t] z[ t] - (f X) z[ t] A 2, Derivative[1][ z] [ t] = = (2 X D4+ D5 Z - c X 2 - e X Z) x[ t] + (D4-2 c X - 
e Z) x[ t] A 2+ (D5+2 X d - e X + f Z) x[ t] z[ t] + (D5 X - p + d X A 2+ f X Z) z[ t] , x[ 0] = =1 .' * A -6, z[ 0] = =1 .' * A -6} , { x, z} , { t, 
0, 30000} , MaxSteps ->50000] 

g001 =Plot[ { X + x[ t] } /. %, { t, 0, 200} , PlotRange -> All, 

PlotStyle -> { RGBColor[ 1 , 0, 0] } , PlotRange -> ALL] 
g003 = Plot[ { Z + z[ t] }/.%%,{ t, 0, 200} , PlotRange -> All, 

PlotStyle -> { RGBColor[ 0, 0, 1] } , PlotRange -> All] 
g004 = ParametricPlot[ Evaluate} { X + x[ t] , Z + z[ t] } / . %%%] , { t, 0, 2000} , 

PlotRange -> All, AxesLabel -> { "X", "Z"} ] 
Show} g001, g003, AxesLabel -> ( "t", "X, Z"} ] 
doi:1 0.1 371 /journal.pone.01 0291 1 .g002 



p«0.8=/> Cj (4.1) 

the concentration of the signaling molecules continuously oscil- 
lates, showing nearly equivalent frequency and amplitude at 
oscillatory individual peaks, except for at the initial phase. 
Subsequently, we evaluated the mean amplitude and frequency 
of the fluctuation oscillation. From t = 50, around when the 
oscillation initiates, to t — 1000, the mean amplitude was calculated 
as the division of the sum of each size of the peak by the number of 
peaks. In addition, the frequency was estimated by dividing the 
number of the peak by At = 950. 

Notably, the present simulation shows that both frequency and 
amplitude are nearly proportional to the logarithm of z=p—p c 
(Figure 3). Namely, 



<f> oc log s (4.2) 

These formulae imply that the outside alteration is transformed 
inside into the information of cell signaling. On the basis of the 
above simulation results, the present signaling model system is 
characterized by: 

7 = log 2 e 
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Figure 4. Scheme of the transformation of outside information into the intracellular signaling oscillation. The scheme illustrates the 
bifurcation of the fluctuation with respect to the increase in p. 
doi:1 0.1 371 /journal.pone.01 0291 1.g004 

Evaluation of the stability of the model around the 
equilibrium state 

For mathematical analysis of stability around the critical point, 
Eq. (3.9) (4.10) was formulated. When p is equivalent to p c = 0.8, 
the matrix for (x, z) is given by Lc in (3.10): 



Further u is given by: 



u = h(e,v) 



(6.5) 



(5.1) 

-{R(Pi -aX e ) + 2D 4 X e +D 5 Z e } (p c -bX e - D 5 X e -dX 2 -fX e Z e ) 
(2D 4 X e + D 5 Z e - cX 2 - eX e Z e ) (D 5 X C + dX 2 +fX e Z e -p c ) 



Using the eigenvectors of L, (//, l 2 ), coordinate transformation is 
performed as follows: 





=[h h] 


u 






z 




_ V 



u 


=[h h]- 1 


X 






_ V _ 




Z 



(6.2) 



(6.3) 



The above parameters are subsequently set as: 

u =/„(«, v) 
v=f v (u,v) 



(6.4) 



= a\ v 2 + a 2 ve + a^e 1 + a 4 v 3 + a 5 v 2 e + a 6 ve 2 + a 7 e ! + 0(4) 

Therefore, using 

dh(e,v) . dh(e,v) . ,„ . . , , 

u= dv vH ^^e = {2aiv + a 2 e)j u (u,v), (6.6) 

we can then obtain using (6.5) and (6.6): 



(2aw + a 2 e)f u (u,v) 

= a 1 v 2 + a 2 ve + a3E 2 + a 4 v 3 + a 5 v 2 e + a 6 ve 2 + a 7 e 3 + 0(4) 



(6.7) 



Solving the above (6.7), the coefficients of a; in (6.5) are given (see 
Appendix SI: mathematica cord [#55-70]). By substitution of u 
that is given by v and E into/^ (u,v) in (6.4), we can obtain the 
kinetic stability equation of fluctuation v using coefficients n, (i = 1, 
2, 3, 4, 5, 6) as follows (See Appendix SI, Out [#73] in the 
mathematica cord): 

v = mv 2 + n 2 EV + n3E 2 v + n4V 2 + n S EV 2 +n6V 3 + 0(4) (6.8) 



To evaluate the amplitude of fluctuation, setting the right hand 
equal to zero, 
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m v 2 + n 2 ve + « 3 e 2 v + n A v 2 + n 5 e v 2 + n 6 v 3 + 0(4) = 0 (6.9) 

As a result, when «<7.0xl0 4 , Eq. (6.9) has two real number 
solutions of li other than zero, indicating the bi-stability of the 
fluctuation v; when ££7.0 x 10 4 (6.9), it has only the zero solution 
of the fluctuation (See Appendix SI; mathematica cord Out 
[7^74]). Therefore, a bifurcation of the fluctuation with respect to 
the value of e is predicted (Figure 4). Because v is simply given by 
the linear equation (6.3), this direcdy demonstrates the amplitude 
bifurcation of x and z. 

Discussion 

Here, we presented a model of cell signaling systems and 
performed mathematical analysis on the model in addition to 
numerical simulations. An increase in the supply of the cofactor 
near the critical concentration induces a 'phase transition' of the 
system, indicating that the model system has the ability to 
transform information on the concentration change of a cofactor 
outside the system into inside information, i.e., the amplitude or 
frequency of the concentration oscillation of the signaling molecule. 
This term, 'phase transition', is a metaphor implying that the model 
system nearly discontinuously acquires the ability to dynamically 
transform outside information into inside information. 

The introduced non-linear kinetic equations include only two 
independent parameters, the active or inactive cofactor binding 
protein. The observed oscillation of signaling molecules in the 
simulation is not a chaotic behavior that requires more than two 
parameters [25]. However, the fluctuation of the signaling 
molecule concentration shows chaos-like oscillatory behavior. In 
fact, neither the amplitude nor frequency of every oscillation is 
precisely constant for a lengthy period as shown in the trajectory. 
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We will report mathematical validation of the result elsewhere. We 
will report mathematical validation of the result elsewhere. 

The simulation allowed us to define the formula (5.3) mentioned 
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